A. Five-step protocol in one dimension.
p_grid <-seq(from=0,to=1,length.out=100)
prior <-rep(1,100)/sum(rep(1,100))
sum(prior)
## [1] 1
likelihood <- dbinom(6,size=9,prob=p_grid)
sum(likelihood)
## [1] 9.900001
unstd.posterior <-likelihood*prior
posterior <-unstd.posterior/sum(unstd.posterior)
sum(posterior)
## [1] 1
plot( p_grid, posterior,
type="b",
xlab="Probability of water",
ylab="Posterior probability")
mtext( "100 points")
p_grid <-seq(from=0,to=1,length.out=100)
prior <-rep(1,100)/sum(rep(1,100))
sum(prior)
## [1] 1
likelihood <-dbinom(6,size=9,prob=p_grid)
sum(likelihood)
## [1] 9.900001
unstd.posterior <-likelihood*prior
posterior <-unstd.posterior/sum(unstd.posterior)
sum(posterior)
## [1] 1
plot( p_grid, posterior,
type="b",
xlab="Probability of water",
ylab="Posterior probability")
mtext( "100 points")
sprior <- ifelse(p_grid < 0.5, 0, 1)/sum(ifelse(p_grid < 0.5, 0, 1))
sum(sprior)
## [1] 1
slikelihood <-dbinom(6,size=9,prob=p_grid)
unstd.sposterior <-slikelihood*sprior
sposterior <-unstd.sposterior/sum(unstd.sposterior)
sum(sposterior)
## [1] 1
plot( p_grid, sposterior,
type="b",
xlab="Probability of water",
ylab="Posterior probability")
points(p_grid, sprior, col='magenta')
points(p_grid, posterior, col='cyan')
mtext( "Step prior, 100 grid points")
legend(x = "topright",
c("sposterior", "sprior", "posterior"),
cex=.8,
col=c("black","magenta","cyan"),
lwd = 2,
bty = "n"
)
bimodal
Pierre Simon LaPlace
Eassy on Probability
There is a newer translation by Andrew I. Dale in 1998 and 2012 published by Springer. It is suppose to be easier to read than the 1902 version and there is commentary on the mathematical notation.
Andrew Dale also wrote
Dale, A.I. (2012) A history of inverse probability: From Thomas Bayes to Karl Pearson Springer Science & Business Media.
and
Dale, A.I. (2003) Most Honourable Remembrance: The Life and Work of Thomas Bayes. Springer.
McGrayne’s book
library(rethinking)
## Loading required package: rstan
## Loading required package: StanHeaders
## Loading required package: ggplot2
## rstan (Version 2.21.5, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## Loading required package: cmdstanr
## This is cmdstanr version 0.5.2
## - CmdStanR documentation and vignettes: mc-stan.org/cmdstanr
## - CmdStan path: /Users/blaine/.cmdstan/cmdstan-2.30.1
## - CmdStan version: 2.30.1
## Loading required package: parallel
## rethinking (Version 2.21)
##
## Attaching package: 'rethinking'
## The following object is masked from 'package:rstan':
##
## stan
## The following object is masked from 'package:stats':
##
## rstudent
globe.qa <- quap(
alist(
W ~ dbinom(W+L, p), # binomial likelihood
p ~ dunif(0,1) # uniform prior
),
data=list(W=6, L=3) )
# display summary of quadratic approximation
precis(globe.qa)
# analytical calculation
W <-6
L <-3
curve( dbeta(x,W+1,L+1),from=0,to=1)
# quadraticapproximation
curve( dnorm(x,0.67,0.16),lty=2,add=TRUE)
globe.qa <- quap(
alist(
W ~ dbinom(W+L, p), # binomial likelihood
p ~ dunif(0,1) # uniform prior
),
data=list(W=12, L=6) )
# display summary of quadratic approximation
precis(globe.qa)
# analytical calculation
W <-12
L <-6
curve( dbeta(x,W+1,L+1),from=0,to=1)
# quadraticapproximation
curve( dnorm(x,0.67,0.11),lty=2,add=TRUE)
globe.qa <- quap(
alist(
W ~ dbinom(W+L, p), # binomial likelihood
p ~ dunif(0,1) # uniform prior
),
data=list(W=24, L=12) )
# display summary of quadratic approximation
precis(globe.qa)
Plug the new sd into the quadratic approximation formula. We finally get more reasonable approximation.
# analytical calculation
W <-24
L <-12
curve( dbeta(x,W+1,L+1),from=0,to=1)
# quadraticapproximation
curve( dnorm(x,0.67,0.08),lty=2,add=TRUE)
globe.qa <- quap(
alist(
W ~ dbinom(W+L, p), # binomial likelihood
p ~ dunif(0,1) # uniform prior
),
data=list(W=48, L=24) )
# display summary of quadratic approximation
precis(globe.qa)
# analytical calculation
W <-48
L <-24
curve( dbeta(x,W+1,L+1),from=0,to=1)
# quadraticapproximation
curve( dnorm(x,0.67,0.06),lty=2,add=TRUE)
Multi-level or hierarchical models are common and have large numbers of parameters. Grid approximation and Quadratic approximation are often not adquate for large models.
MCMC merely draws samples from the posterior rather than computing the posterior directly. The frequencies of the parameter values of these samples is proportional to their posterior plausability. The histogram of these samples gives a picture of the posterior. So with MCMC, we work with samples rather then with an estimate of the posterior distribution.
Source: Sobol, I. M. (1973). Numerical Monte Carlo Methods. Nauka, Moscow.
The ratio of the number of points inside the boundary and outside give the proportion of the square occupied by the irregular area.
#### The bounding area could
be another curve
Source: Robert, C. P., Casella, G., and Casella, G. (2010). Introducing Monte Carlo methods with R. New York: Springer.
Determine area under aribtrary curve
\[\begin{equation} P(\Delta E)=e^{-\Delta E / k T} \end{equation}\]
Make a move (e.g., select a sample, change a torsion angle, translate the ligand, rotate the ligand)
Calculate the energy (\(E\))
Compare E to prior \(E^{\circ}\)
Source: Metropolis, Nicholas; Rosenbluth, Arianna W.; Rosenbluth, Marshall N.; Teller, Augusta H.; and Teller, Edward (1953) Equation of state Calculations by Fast Computing Machines. The Journal of Chemical Physics. 21(6):1087.
n_samples <-1000
p <-rep(NA,n_samples)
p[1] <- 0.5
W <-6
L <-3
for (i in 2:n_samples){
p_new <-rnorm(1,p[i-1],0.1)
if (p_new<0) p_new <- abs(p_new)
if (p_new>1) p_new <- 2 - p_new
q0 <- dbinom(W,W+L,p[i-1])
q1 <- dbinom(W,W+L,p_new)
p[i] <- ifelse(runif(1)<q1/q0,p_new,p[i-1])
}
dens( p, xlim=c(0,1))
curve( dbeta(x,W+1,L+1),lty=2,add=TRUE)
n_samples <-10000
p <-rep(NA,n_samples)
p[1] <- 0.5
W <-6
L <-3
for (i in 2:n_samples){
p_new <-rnorm(1,p[i-1],0.1)
if (p_new<0) p_new <- abs(p_new)
if (p_new>1) p_new <- 2 - p_new
q0 <- dbinom(W,W+L,p[i-1])
q1 <- dbinom(W,W+L,p_new)
p[i] <- ifelse(runif(1)<q1/q0,p_new,p[i-1])
}
dens( p, xlim=c(0,1))
curve( dbeta(x,W+1,L+1),lty=2,add=TRUE)
n_samples <-100000
p <-rep(NA,n_samples)
p[1] <- 0.5
W <-6
L <-3
for (i in 2:n_samples){
p_new <-rnorm(1,p[i-1],0.1)
if (p_new<0) p_new <- abs(p_new)
if (p_new>1) p_new <- 2 - p_new
q0 <- dbinom(W,W+L,p[i-1])
q1 <- dbinom(W,W+L,p_new)
p[i] <- ifelse(runif(1)<q1/q0,p_new,p[i-1])
}
dens( p, xlim=c(0,1))
curve( dbeta(x,W+1,L+1),lty=2,add=TRUE)
n_samples <-1000000
p <-rep(NA,n_samples)
p[1] <- 0.5
W <-6
L <-3
for (i in 2:n_samples){
p_new <-rnorm(1,p[i-1],0.1)
if (p_new<0) p_new <- abs(p_new)
if (p_new>1) p_new <- 2 - p_new
q0 <- dbinom(W,W+L,p[i-1])
q1 <- dbinom(W,W+L,p_new)
p[i] <- ifelse(runif(1)<q1/q0,p_new,p[i-1])
}
dens( p, xlim=c(0,1))
curve( dbeta(x,W+1,L+1),lty=2,add=TRUE)